M2.855 · Modelos avanzados de minería de datos · PEC2
2021-1 · Máster universitario en Ciencia de datos (Data science)
Estudios de Informática, Multimedia y Telecomunicación
A lo largo de esta práctica veremos como aplicar distintas técnicas no supervisadas así como algunas de sus aplicaciones reales:
Para ello vamos a necesitar las siguientes librerías:
import random
import umap
import numpy as np
import pandas as pd
import sklearn
from sklearn import cluster # Algoritmos de clustering.
from sklearn import datasets # Crear datasets.
from sklearn import decomposition # Algoritmos de reduccion de dimensionalidad.
# Visualizacion.
import matplotlib
import matplotlib.pyplot as plt
%matplotlib inline
Este ejercicio trata de explorar distintas técnicas de agrupamiento ajustándolas a distintos conjuntos de datos.
El objetivo es doble: entender la influencia de los parámetros en su comportamiento, y conocer sus limitaciones en la búsqueda de estructuras de datos.
X_blobs, y_blobs = datasets.make_blobs(n_samples=1000, n_features=2, centers=4, cluster_std=1.6, random_state=42)
X_moons, y_moons = datasets.make_moons(n_samples=1000, noise=.07, random_state=42)
X_circles, y_circles = datasets.make_circles(n_samples=1000, factor=.5, noise=.05, random_state=42)
Cada dataset tiene 2 variables: una variable X que contiene 2 features (columnas) y tantas filas como muestras. Y una variable y que alberga las etiquetas que identifican cada cluster.
A lo largo del ejercicio no se usará la variable y (sólo con el objetivo de visualizar). El objetivo es a través de los distintos modelos de clustering conseguir encontrar las estructuras descritas por las variables y.
fig, axis = plt.subplots(2, 3, figsize=(13, 8))
for i, (X, y, ax, name) in enumerate(zip([X_blobs, X_moons, X_circles] * 2,
[None] * 3 + [y_blobs, y_moons, y_circles],
axis.reshape(-1),
['Blob', 'Moons', 'Circles'] * 2)):
ax.set_title('Dataset: {}, '.format(name) + ('lo que analizarás' if i < 3 else 'los grupos a encontrar'))
ax.scatter(X[:,0], X[:,1], s=15, c=y, alpha=.3, cmap='jet')
ax.axis('equal')
plt.tight_layout()
En este apartado se pide probar el algoritmo k-means sobre los tres datasets presentados anteriormente ajustando con los parámetros adecuados y analizar sus resultados.
X, y = X_blobs, y_blobs
Para estimar el número de clusters a detectar por k-means. Una técnica para estimar $k$ es, como se explica en la teoría:
Los criterios anteriores (minimización de distancias intra grupo o maximización de distancias inter grupo) pueden usarse para establecer un valor adecuado para el parámetro k. Valores k para los que ya no se consiguen mejoras significativas en la homogeneidad interna de los segmentos o la heterogeneidad entre segmentos distintos, deberían descartarse.
Lo que popularmente se conocer como regla del codo.
Primero es necesario calcular la suma de los errores cuadráticos (SSE) que consiste en la suma de todos los errores (distancia de cada punto a su centroide asignado) al cuadrado.
$$SSE = \sum_{i=1}^{K} \sum_{x \in C_i} euclidean(x, c_i)^2$$Donde $K$ es el número de clusters a buscar por k-means, $x \in C_i$ son los puntos que pertenecen a i-ésimo cluster, $c_i$ es el centroide del cluster $C_i$ (al que pertenece el punto $x$), y $euclidean$ es la distancia euclídea.
Este procedimiento realizado para cada posible valor $k$, resulta en una función monótona decreciente, donde el eje $x$ representa los distintos valores de $k$, y el eje $y$ el $SSE$. Intuitivamente se podrá observar un significativo descenso del error, que indicará el valor idóneo de $k$.
Se pide realizar la representación gráfica de la regla del codo junto a su interpretación, utilizando la librería matplotlib y la implementación en scikit-learn de k-means.
#############################################
# SOLUCION #
#############################################
from typing import List
def plot_elbow_rule(datasets: List[np.array],
min_k: int = 1,
max_k: int = 10) -> None:
"""Computes and shows SSE for each value of k in k-means.
This is a generic solution that receives a list containing datasets,
and shows a plot for each datset.
Args:
datasets: list
A list of np.array with shape (points, dimensions).
min_k: int
The lowest value of k to try.
max_k: int
The highest value of k to try.
"""
fig, ax = plt.subplots(1, len(datasets), figsize=(8, 4))
for i, a in zip(range(len(datasets)), ax if type(ax) == np.ndarray else (ax,)):
x = range(min_k, max_k + 1)
y = []
X = datasets[i]
for nc in x:
kmeans = cluster.KMeans(n_clusters=nc)
dists = kmeans.fit_transform(X)
y.append(np.sum(np.min(dists, axis=1) ** 2))
a.plot(x, y)
a.set_xlabel('k (# de clusters)')
a.set_ylabel('SSE')
a.set_title('Dataset {}'.format(i))
a.set_xticks(range(min_k, max_k + 1))
a.grid()
plt.tight_layout()
dists = plot_elbow_rule([X])
#############################################
# SOLUCION #
#############################################
fig, ax = plt.subplots(1, 1, figsize=(5, 5))
kmeans = cluster.KMeans(n_clusters=4, algorithm='full')
ax.scatter(X[:,0], X[:,1], c=np.argmin(kmeans.fit_transform(X), axis=1), alpha=.3, cmap='jet')
ax.scatter(kmeans.cluster_centers_[:,0], kmeans.cluster_centers_[:,1],
s=350, c='k', marker='X', linewidths=2, edgecolors='w')
ax.axis('equal')
plt.tight_layout()
X, y = X_moons, y_moons
#############################################
# SOLUCION #
#############################################
dists = plot_elbow_rule([X])
#############################################
# SOLUCION #
#############################################
fig, ax = plt.subplots(1, 1, figsize=(5, 5))
kmeans = cluster.KMeans(n_clusters=2, algorithm='full')
ax.scatter(X[:,0], X[:,1], c=np.argmin(kmeans.fit_transform(X), axis=1), alpha=.3, cmap='jet')
ax.scatter(kmeans.cluster_centers_[:,0], kmeans.cluster_centers_[:,1],
s=350, c='k', marker='X', linewidths=2, edgecolors='w')
ax.axis('equal')
plt.tight_layout()
X, y = X_circles, y_circles
#############################################
# SOLUCION #
#############################################
dists = plot_elbow_rule([X])
#############################################
# SOLUCION #
#############################################
fig, ax = plt.subplots(1, 1, figsize=(5, 5))
kmeans = cluster.KMeans(n_clusters=2, algorithm='full')
ax.scatter(X[:,0], X[:,1], c=np.argmin(kmeans.fit_transform(X), axis=1), alpha=.3, cmap='jet')
ax.scatter(kmeans.cluster_centers_[:,0], kmeans.cluster_centers_[:,1],
s=350, c='k', marker='X', linewidths=2, edgecolors='w')
ax.axis('equal')
plt.tight_layout()
X, y = X_blobs, y_blobs
#############################################
# SOLUCION #
#############################################
model = cluster.DBSCAN(eps=2.2, min_samples=100, n_jobs=-1)
clusters = model.fit(X)
fig, ax = plt.subplots(1, 1, figsize=(5, 5))
ax.scatter(X[:,0], X[:,1], c=clusters.labels_, alpha=.3, cmap='jet')
ax.axis('equal')
plt.tight_layout()
X, y = X_moons, y_moons
#############################################
# SOLUCION #
#############################################
model = cluster.DBSCAN(eps=0.15, min_samples=10, n_jobs=-1)
clusters = model.fit(X)
fig, ax = plt.subplots(1, 1, figsize=(5, 5))
ax.scatter(X[:,0], X[:,1], c=clusters.labels_, alpha=.3, cmap='jet')
ax.axis('equal')
plt.tight_layout()
X, y = X_circles, y_circles
#############################################
# SOLUCION #
#############################################
model = cluster.DBSCAN(eps=0.12, min_samples=5, n_jobs=-1)
clusters = model.fit(X)
fig, ax = plt.subplots(1, 1, figsize=(5, 5))
ax.scatter(X[:,0], X[:,1], c=clusters.labels_, alpha=.3, cmap='jet')
ax.axis('equal')
plt.tight_layout()
En este apartado se pide visualizar mediante un dendrograma la construcción progresiva de los grupos mediante un algoritmo jerárquico aglomerativo (estrategia bottom-up). Con ello se pretende encontrar un método gráfico para entender el comportamiento del algoritmo y encontrar los clusters deseados en cada dataset.
X, y = X_blobs, y_blobs
#############################################
# SOLUCION #
#############################################
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
fig, ax = plt.subplots(1, 2, figsize=(11, 5))
lin = linkage(pdist(X), 'complete')
dendrogram(lin, ax=ax[0])
ax[0].set_title('Dataset {} con linkage "{}"'.format(X, 'complete'))
cluster_idx = fcluster(lin, t=15, criterion='distance')
ax[1].scatter(X[:, 0], X[:, 1], c=cluster_idx, alpha=.3, cmap='jet')
ax[1].axis('equal')
plt.tight_layout()
#############################################
# SOLUCION #
#############################################
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
fig, ax = plt.subplots(1, 2, figsize=(11, 5))
lin = linkage(pdist(X), 'ward')
dendrogram(lin, ax=ax[0])
ax[0].set_title('Dataset {} con linkage "{}"'.format(X, 'ward'))
cluster_idx = fcluster(lin, t=50, criterion='distance')
ax[1].scatter(X[:, 0], X[:, 1], c=cluster_idx, alpha=.3, cmap='jet')
ax[1].axis('equal')
plt.tight_layout()
X, y = X_moons, y_moons
#############################################
# SOLUCION #
#############################################
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
fig, ax = plt.subplots(1, 2, figsize=(11, 5))
lin = linkage(pdist(X), 'single')
dendrogram(lin, ax=ax[0])
ax[0].set_title('Dataset {} con linkage "{}"'.format(X, 'single'))
cluster_idx = fcluster(lin, t=.2, criterion='distance')
ax[1].scatter(X[:, 0], X[:, 1], c=cluster_idx, alpha=.3, cmap='jet')
ax[1].axis('equal')
plt.tight_layout()
X, y = X_circles, y_circles
#############################################
# SOLUCION #
#############################################
from scipy.spatial.distance import pdist
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
fig, ax = plt.subplots(1, 2, figsize=(11, 5))
lin = linkage(pdist(X), 'single')
dendrogram(lin, ax=ax[0])
ax[0].set_title('Dataset {} con linkage "{}"'.format(X, 'single'))
cluster_idx = fcluster(lin, t=.2, criterion='distance')
ax[1].scatter(X[:, 0], X[:, 1], c=cluster_idx, alpha=.3, cmap='jet')
ax[1].axis('equal')
plt.tight_layout()
Es posible aplicar una amplia variedad de algoritmos para la reducción de dimensionalidad. Para ello se empleará el dataset MNIST compuesto de miles de dígitos manuscritos del 0 al 9. Donde cada imagen se compone de 784 píxeles (imágenes de 28 x 28), por lo que se parte de un número alto de dimensiones.
X, y = sklearn.datasets.fetch_openml('mnist_784', version=1, return_X_y=True, as_frame=False)
Por lo que cada muestra (las 70k filas del dataset) se componen de 784 dimensiones:
X.shape
(70000, 784)
fig, axis = plt.subplots(1, 10, figsize=(12, 6))
for i, ax in enumerate(axis):
ax.imshow(X[y == str(i)][0].reshape(28, 28), cmap='gray_r')
ax.set_title(str(i))
ax.axis('off')
plt.tight_layout()
Si cada algoritmo obtiene resultados distintos a la hora de reducir la dimensionalidad, ¿qué representación es más fiel a la distribución original?
Antes de reducir las 784 dimensiones originales de cada muestra a 2 para poder visualizarlas en 2 dimensiones, es muy útil conocer, o al menos intuir, la estructura en alta dimensionalidad de los datos.
Para ello se puede hacer uso del dendrograma como heurística para conocer la disposición original de los datos y comprobar si la proyección es similar a lo mostrado por el dendrograma.
#############################################
# SOLUCION #
#############################################
model = decomposition.PCA(n_components=2, random_state=42)
model.fit(X)
X_t = model.transform(X)
fig, ax = plt.subplots(1, 1, figsize=(10, 8))
for i in range(10):
ax.scatter(X_t[y == str(i),0], X_t[y == str(i),1], s=3, label=str(i), alpha=.2)
plt.legend()
plt.tight_layout()
En la gráfica anterior cada punto representa una muestra en 2 dimensiones. Con PCA es posible invertir la transformación para que, a partir de cada punto 2d, se obtenga de nuevo (aproximadamente) la imagen original (784 dimensiones).
Por lo que es posible "generar" nuevas imágenes eligiendo puntos aleatoriamente del plano 2d, y pedirle al modelo PCA aprendido que invierta la transformación para obtener las "teóricas" imágenes que habrían sido proyectadas a esos puntos del espacio proyectado.
#############################################
# SOLUCION #
#############################################
min_ = X_t.min(axis=0)
max_ = X_t.max(axis=0)
range_0 = np.linspace(min_[0], max_[0], num=10)
range_1 = np.linspace(min_[1], max_[1], num=10)
Con las dos secuencias de 10 (una por cada dimensión del plano de proyección) valores es posible combinar los puntos de ambas secuencias para generar 100 combinaciones (puntos 2d) que teselan el plano sobre el que PCA ha proyectado las muestras.
#############################################
# SOLUCION #
#############################################
fig, ax = plt.subplots(10, 10, figsize=(10, 10))
for i, c0 in enumerate(range_0):
for j, c1 in enumerate(range_1):
ax[i, j].imshow(model.inverse_transform([c0, c1]).reshape(28, 28), cmap='gray_r')
ax[i, j].set_yticks([])
ax[i, j].set_xticks([])
plt.tight_layout()
#############################################
# SOLUCION #
#############################################
model = umap.UMAP(n_components=2, random_state=42)
model.fit(X)
X_t = model.transform(X)
fig, ax = plt.subplots(1, 1, figsize=(10, 8))
for i in range(10):
ax.scatter(X_t[y == str(i),0], X_t[y == str(i),1], s=3, label=str(i), alpha=.2)
plt.legend()
plt.tight_layout()
OMP: Info #276: omp_set_nested routine deprecated, please use omp_set_max_active_levels instead.
#############################################
# SOLUCION #
#############################################
min_ = X_t.min(axis=0)
max_ = X_t.max(axis=0)
range_0 = np.linspace(min_[0], max_[0], num=12)[1:-1]
range_1 = np.linspace(min_[1], max_[1], num=12)[1:-1]
l = []
for i, c0 in enumerate(range_0):
for j, c1 in enumerate(range_1):
l.append([c0, c1])
out = model.inverse_transform(l)
fig, ax = plt.subplots(10, 10, figsize=(12, 10))
c = 0
for i in range(10):
for j in range(10):
ax[i, j].imshow(out[c].reshape(28, 28), cmap='gray_r')
ax[i, j].set_yticks([])
ax[i, j].set_xticks([])
c += 1
En este ejercicio se busca automatizar la localización de lugares turísticos a través de los metadatos de las fotografías de flickr.
Para ello se provee junto a la PEC el dataset: barcelona.csv. Ya que se pide encontrar los puntos de mayor interés turístico de esta ciudad.
Opcional: si quieres hacerlo para otra región
Pero si quieres hacerlo para otra parte del mundo, puedes descargarte el dataset completo aquí y descomprime para extraer el CSV.
Para seleccionar las coordenadas de la zona de interés puedes usar la opción Export manual de OpenStreetMaps.
Por último, para filtrar los datos que se corresponden a la zona deseada puedes usar el programa AWK mediante la siguiente línea:
awk -F"," 'NR == 1 {print $5","$6} (NR > 1 && $5 > 41.3560 && $5 < 41.4267 && $6 > 2.1300 && $6 < 2.2319) {print $5","$6}' photo_metadata.csv
$5 hace referencia a la latitud, y $6 a la longitud. Sustituye los valores mínimo y máximo para obtener los datos de localización referentes a tu área de interés.
geo_df = pd.read_csv('barcelona.csv', header=0)
geo_df.sample(5)
| latitude | longitude | |
|---|---|---|
| 2382 | 41.386500 | 2.155999 |
| 17594 | 41.383504 | 2.181404 |
| 6288 | 41.418015 | 2.153320 |
| 17535 | 41.379658 | 2.174177 |
| 892 | 41.395288 | 2.161770 |
#############################################
# SOLUCION #
#############################################
fig, ax = plt.subplots(1, 1, figsize=(10, 7))
ax.scatter(geo_df['longitude'], geo_df['latitude'], s=5, alpha=.3)
ax.axis('equal')
plt.tight_layout()
#############################################
# SOLUCION #
#############################################
geo_df_sample = geo_df.sample(frac=.5)
fig, ax = plt.subplots(1, 1, figsize=(10, 7))
ax.scatter(geo_df_sample['longitude'], geo_df_sample['latitude'], s=5, alpha=.3)
ax.axis('equal')
plt.tight_layout()
#############################################
# SOLUCION #
#############################################
X = geo_df_sample.values
model = cluster.DBSCAN(eps=0.0011, min_samples=40, n_jobs=-1)
clusters = model.fit(X)
fig, ax = plt.subplots(1, 1, figsize=(10, 7))
ax.scatter(X[:,1], X[:,0], s=5, c=clusters.labels_, alpha=.3, cmap='tab20')
ax.set_title('{} clusters encontrados'.format(clusters.labels_.max()))
ax.axis('equal')
plt.tight_layout()
#############################################
# SOLUCION #
#############################################
fig, ax = plt.subplots(1, 1, figsize=(10, 7))
ax.set_title('{} clusters encontrados'.format(clusters.labels_.max()))
ax.scatter(X[:,1][clusters.labels_ >= 0], X[:,0][clusters.labels_ >= 0], s=5,
c=clusters.labels_[clusters.labels_ >= 0], alpha=.3, cmap='tab20')
ax.axis('equal')
plt.tight_layout()
import smopy
map = smopy.Map((41.3560, 2.1300, 41.4267, 2.2319), z=15)
ax = map.show_mpl(figsize=(20, 15))
for i in range(clusters.labels_.max() + 1):
x, y = map.to_pixels(X[:,0][clusters.labels_ == i].mean(), X[:,1][clusters.labels_ == i].mean())
ax.plot(x, y, 'or', ms=10, mew=2);
Lowered zoom level to keep map size reasonable. (z = 13)